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Summary. In this paper we consider two sources of enhancement for the meshfree 
Lagrangian particle method smoothed particle hydrodynamics (SPH) by improving 
the accuracy of the particle approximation. Namely, we will consider shape functions 
constructed using: moving least-squares approximation (MLS); radial basis functions 
(RBF). Using MLS approximation is appealing because polynomial consistency of 
the particle approximation can be enforced. RBFs further appeal as they allow one 
to dispense with the smoothing-length - the parameter in the SPH method which 
governs the number of particles within the support of the shape function. Currently, 
only ad hoc methods for choosing the smoothing-length exist. We ensure that any 
enhancement retains the conservative and meshfree nature of SPH. In doing so, 
we derive a new set of variationally-consistent hydrodynamic equations. Finally, we 
demonstrate the performance of the new equations on the Sod shock tube problem. 



1 Introduction 

Smoothed particle hydrodynamics (SPH) is a meshfree Lagrangian particle 
method primarily used for solving problems in solid and fluid mechanics 
(see [10] for a recent comprehensive review). Some of the attractive character- 
istics that SPH possesses include: the ability to handle problems with large 
deformation, free surfaces and complex geometries; truly meshfree nature (no 
background mesh required) ; exact conservation of momenta and total energy. 
On the other hand, SPH suffers from several drawbacks: an instability in 
tension; difficulty in enforcing essential boundary conditions; fundamentally 
based on inaccurate kernel approximation techniques. This paper addresses 
the last of these deficiencies by suggesting improved particle approximation 
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procedures. Previous contributions in this direction (reviewed in [2]) have fo- 
cused on corrections of the existing SPH particle approximation (or its deriva- 
tives) by enforcing polynomial consistency. As a consequence, the conservation 
of relevant physical quantities by the discrete equations is usually lost. 

The outline of the paper is as follows. In the next section we review how 
SPH equations for the non-dissipative motion of a fluid can be derived. In 
essence this amounts to a discretization of the Euler equations: 

dp „ dv 1_,t, de P_ 

■£ = -pV-v, — — VP, — = V-w, 1 

dt dt p dt p 

where ^ is the total derivative, p, v, e and P are the density, velocity, ther- 
mal energy per unit mass and pressure, respectively. The derivation is such 
that important conservation properties are satisfied by the discrete equations. 
Within the same section we derive a new set of variationally-consistent hy- 
drodynamic equations based on improved particle approximation. In Sect. 3 
we construct specific examples - based on moving least-squares approxima- 
tion and radial basis functions - to complete the newly derived equations. 
The paper finishes with Sect. 4 where we demonstrate the performance of the 
new methods on the Sod shock tube problem [12] and make some concluding 
remarks. 

To close this section, we briefly review the SPH particle approximation 
technique on which the SPH method is fundamentally based and which 
we purport to be requiring improvement. From a set of scattered particles 
{xi, . . . , xjy} C fft d , SPH particle approximation is achieved using 

N 

Sf(x) = J2f(^W(\x- Xj \,h), (2) 

3 = 1 Pj 

where mj and pj denotes the mass and density of the jth particle, respectively. 
The function W is a normalised kernel function which approximates the 6- 
distribution as the smoothing-length, h, tends to zero. The function ^-W (\x — 
xj \,h) is called an SPH shape function and the most popular choice for W is 
a compactly supported cubic spline kernel with support 2h. The parameter 
h governs the extent to which contributions from neighbouring particles are 
allowed to smooth the approximation to the underlying function /. Allowing a 
spatiotemporally varying smoothing-length increases the accuracy of an SPH 
simulation considerably. There are a selection of ad hoc techniques available 
to accomplish this, although often terms arising from the variation in h are 
neglected in the SPH method. The approximating power of the SPH particle 
approximation is perceived to be poor. The SPH shape functions fail to provide 
a partition of unity so that even the constant function is not represented 
exactly. There is currently no approximation theory available for SPH particle 
approximation when the particles are in general positions. The result of a 
shock tube simulation using the SPH equations derived in Sect. 2 is shown in 
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Fig. 1. Shock tube simulation (t=0.2) using SPH. 

Fig. 1 (see Sect. 4 for the precise details of the simulation). The difficulty that 
SPH has at the contact discontinuity (x ~ 0.2) and the head of the rarefaction 
wave (x « —0.25) is attributed to a combination of the approximation (2) and 
the variable smoothing-lcngth not being sclf-consistently incorporated. 

2 Variationally-Consistent Hydrodynamic Equations 

It is well known (see [10] and the references cited therein) that the most 
common SPH equations for the non-dissipative motion of a fluid can be derived 
using the Lagrangian for hydrodynamics and a variational principle. In this 
section we review this procedure for a particular formulation of SPH before 
deriving a general set of variationally-consistent hydrodynamic equations. 

The aforementioned Lagrangian is a particular functional of the dynamical 
coordinates: L(x, v) = J p(v 2 /2 — e) dx, where x is the position, v is the 
velocity, p is the density, e is the thermal energy per unit mass and the integral 
is over the volume being discretized. Given N particles {x\ . . . , x^} C H d , 
the SPH discretization of the Lagrangian, also denoted by L, is given by 

N v 2 

L = ^l m o{-f - e i)> ( 3 ) 
i=i 

where rrij has replaced PjVj to denote particle mass (assumed to be constant), 
and Vj is a volume associated with each particle. Self-evidently, the notation 
fj is used to denote the function / evaluated at the jth particle. 

The Eulcr-Lagrangc equations give rise to SPH equations of motion pro- 
vided each quantity in (3) can be written directly as a function of the particle 
coordinates. By setting / = p in (2) and evaluating at xj, we can obtain 
an expression for pj directly as a function of the particle coordinates. There- 
fore, because we assume that ej = ej(pj), the Euler-Lagrange equations are 
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amenable. Furthermore, in using this approach, conservation of momenta and 
total energy are guaranteed via Noether's symmetry theorem. However, when 
we consider improved particle approximation, the corresponding expression 
for density depends on the particle coordinates in an implicit manner, so that 
the Euler-Lagrange equations are not directly amenable. To circumvent this 
difficulty, one can use the principle of stationary action directly to obtain SPH 
equations of motion - the action, 



-I 



Ldt, 



being the time integral of L. The principle of stationary action demands that 
the action is invariant with respect to small changes in the particle coordi- 
nates (i.e., SS — 0). The Euler-Lagrange equations are a consequence of this 
variational principle. In [10] it is shown that if an expression for the time rate 
of change of pj is available, then, omitting the detail, this variational principle 
gives rise to SPH equations of motion. 

To obtain an expression for the time rate of change of density we can 
discretize the first equation of (1) using (2) by collocation. By assuming that 
the SPH shape functions form a partition of unity we commit error but are 
able to artificially provide the discretization with invariance to a constant 
shift in velocity (Galilean invariance): 



^ r = -PiY j r ? ± {v j -v i )-V i W{\x i -x j \,h i ), i = l,...,N, (4) 
j=i Pi 

where V, is the gradient with respect to the coordinates of the ith particle. 
The equations of motion that are variationally-consistent with (4) are 

= _]_J2?!li(p i v i w(\x i - Xj \,hi) + PjViWilxi - Xjlhjj), (5) 
at pi pj \ J 

for i = 1, . . . , N, where Pi denotes the pressure of the ith particle (provided 
via a given equation of state). Using the first law of thermodynamics, the 
equation for the rate of change of thermal energy is given by 



dei /', dp, . , T 



As already noted, a beneficial consequence of using the Euler-Lagrange 
equations is that one automatically preserves, in the discrete equations, fun- 
damental conservation properties of the original system (1). Since we have 
not done this, conservation properties are not necessarily guaranteed by our 
discrete equations (4)-(6). However, certain features of the discretization (4) 
give us conservation. Indeed, by virtue of (4) being Galilean invariant, one 
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conserves linear momentum and total energy (assuming perfect time integra- 
tion). Remember that Galilean invariance was installed under the erroneous 
assumption that the SPH shape functions provide a partition of unity. Angu- 
lar momentum is also explicitly conserved by this formulation due to W being 
symmetric. 

Now, we propose to enhance SPH by improving the particle approxima- 
tion (2). Suppose we have constructed shape functions <f>j that provide at least 
a partition of unity. With these shape functions we form a quasi- interpolant: 

N 

Sf = 52f(xj)4>j, (7) 

which we implicitly assume provides superior approximation quality than that 
provided by (2) . We defer particular choices for <frj until the next section. The 
discretization of the continuity equation now reads 

dp N 

-j^ = -Pi^2(vj -Vi) -V<f>j{xi), i = l,...,N, (8) 

where, this time, we have supplied genuine Galilean invariance, without com- 
mitting an error, using the partition of unity property of <pj. As before, the 
principle of stationary action provides the equations of motion and conser- 
vation properties of the resultant equations reflect properties present in the 
discrete continuity equation (8). 

To obtain (3), two assumptions were made. Firstly, the SPH shape func- 
tions were assumed to form a sufficiently good partition of unity. Secondly, 
it was assumed that the kernel approximation J fW(\ ■ —Xj\,h)dx « f(xj), 
was valid. For our general shape functions the first of these assumptions is 
manifestly true. The analogous assumption we make to replace the second is 
that the error induced by the approximations 

dx « ./• J cf> 3 dx^fjVj, j 1 V. (9) 

is negligible. With the assumption (9), the approximate Lagrangian associated 
with <j>j is identical in form to (3). Neglecting the details once again, which 
can be recovered from [10], the equations of motion variationally-consistent 
with (8) are 

j 1 N 

j=i j 

The equations (6), (8) and (10) constitute a new set of variationally- 
consistent hydrodynamic equations. They give rise to the formulation of SPH 
derived earlier under the transformation (f>j(xi) \— > ^-W(\xi — Xj\,hi). The 
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equations of motion (10) appear in [8] but along side variationally- inconsistent 
companion equations. The authors advocate using a variationally-consistcnt 
set of equation because evidence from the SPH literature (e.g., [3, 9]) suggests 
that not doing so can lead to poor numerical results. 

Linear momentum and total energy are conserved by the new equations, 
and this can be verified immediately using the partition of unity property of 
<j>j. The (f>j will not be symmetric. However, if it is also assumed that the 
shape functions reproduce linear polynomials, namely, Yl x j < t > j(, x ) — x, then 
it is simple to verify that angular momentum is also explicitly conserved. 

3 Moving Least-Squares and Radial Basis Functions 

In this section we construct quasi- interpolants of the form (7). In doing so 
we furnish our newly derived hydrodynamic equations (6), (8) and (10) with 
several examples. 

Moving least-squares (MLS). The preferred construction for MLS shape func- 
tions, the so-called Backus-Gilbert approach [4], seeks a quasi- interpolant of 
the form (7) such that: 

• Sp = p for all polynomials p of some fixed degree; 

• <j)j(x), j — 1, . . . , N, minimise the quadratic form ^ (j)'j(x) 

where mi is a fixed weight function. If w is continuous, compactly supported 
and positive on its support, this quadratic minimisation problem admits a 
unique solution. Assuming / has sufficient smoothness, the order of conver- 
gence of the MLS approximation (7) directly reflects the degree of polynomial 
reproduced [14]. 

The use of MLS approximation in an SPH context has been considered be- 
fore. Indeed, Bclytschko et al. [2] have shown that correcting the SPH particle 
approximation up to linear polynomials is equivalent to an MLS approxima- 
tion with w(\ ■ —Xj\/h) = W(\ ■ —Xj\, h). There is no particular reason to base 
the MLS approximation on an SPH kernel. We find that MLS approximations 
based on Wendland functions [13], which have half the natural support of a 
typical SPH kernel, produce results which are less noisy. Dilts [7, 8] employs 
MLS approximation too. Indeed, in [7], Dilts makes an astute observation that 
addresses an inconsistency that arises due to (9) - we have the equations 

= V i J2( V J - v i) ■ Wi(^) and ^ ~ ^ (J M x ) dx j ■ 

Dilts shows that if hi is evolved according to hi oc V^ d then there is agreement 
between the right-hand sides of these equations when a one-point quadra- 
ture of J <pi Ax is employed. Thus, providing some theoretical justification for 
choosing this particular variable smoothing-length over other possible choices. 
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Radial basis functions (RBFs). To construct an RBF interpolant to an un- 
known function / on xi, . . . , Xn, one produces a function of the form 



where the Xj are found by solving the linear system If(xi) — f(xi), i — 
1,...,N. The radial basis function, ip, is a pre-specified univariate function 
chosen to guarantee the solvability of this system. Depending on the choice of 
ip, a low degree polynomial is sometimes added to (11) to ensure solvability, 
with the additional degrees of freedom taken up in a natural way. This is 
the case with the polyharmonic splines, which are defined, for m > d/2, by 
■0(M) = |a;| 2m_d log \x\ if d is even and VKM) = |a;| 2m_<i otherwise, and a 
polynomial of degree m — 1 is added. The choice m > 2 ensures the RBF 
interpolant reproduces linear polynomials as required for angular momentum 
to be conserved by the equations of motion. As with MLS approximation, 
one has certain strong assurances regarding the quality of the approximation 
induced by the RBF interpolant (e.g. [6] for the case of polyharmonic splines). 

In its present form (11), the RBF interpolant is not directly amenable. 
One possibility is to rewrite the interpolant in cardinal form so that it co- 
incides with (7). This naively constitutes much greater computational effort. 
However, there are several strategies for constructing approximate cardinal 
RBF shape functions (e.g. [5]) and fast evaluation techniques (e.g. [1]) which 
reduce this work significantly. The perception of large computational effort is 
an attributing factor as to why RBFs have not been considered within an SPH 
context previously. Specifically for polyharmonic splines, another possibility 
is to construct shape functions based on discrete m-iterated Laplacians of ip. 
This is sensible because the continuous iterated Laplacian, when applied ip, 
results in the (5-distribution (up to a constant). This is precisely the approach 
we take in Sect. 4 where we employ cubic B-spline shape functions for one of 
our numerical examples. The cubic B-splines are discrete bi-Laplacians of the 
shifts of | • | 3 , and they gladly reproduce linear polynomials. 

In addition to superior approximation properties, using globally supported 
RBF shape functions has a distinct advantage. One has dispensed with the 
smoothing-length entirely. Ducly, issues regarding how to correctly vary and 
self-consistently incorporate the smoothing-length vanish. Instead, a natural 
'support' is generated related to the relative clustering of particles. 

4 Numerical Results 

In this section we demonstrate the performance of the scheme (6), (8) and (10) 
using both MLS and RBF shape functions. The test we have selected has 
become a standard one-dimensional numerical test in compressible fluid flow 
- the Sod shock tube [12]. The problem consists of two regions of ideal gas, 



N 




(11) 
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one with a higher pressure and density than the other, initially at rest and 
separated by a diaphragm. The diaphragm is instantaneously removed and the 
gases allowed to flow resulting in a rarefaction wave, contact discontinuity and 
shock. We set up 450 equal mass particles in [—0.5,0.5]. The gas occupying 
the left-hand and right-hand sides of the domain arc given initial conditions 
(Pl,Pl,v l ) = (1.0,1.0,0.0) and {P R ,p R ,v R ) = (0.1,0.125,0.0), respectively. 
The initial condition is not smoothed. 

With regards to implementation, artificial viscosity is included to prevent 
the development of unphysical oscillations. The form of the artificial viscosity 
mimics that of the most popular SPH artificial viscosity and is applied with a 
switch which reduces the magnitude of the viscosity by a half away from the 
shock. A switch is also used to administer an artificial thermal conductivity 
term, also modelled in SPH. Details of both dissipative terms and their re- 
spective switches can be accessed through [10]. Finally, we integrate, using a 
predictor-corrector method, the equivalent hydrodynamic equations 

dV N 

N 

d Vl 1 ^ de 4 P l dV t 

together with ^ = Vi, to move the particles. To address the consistency issue 
regarding particle volume mentioned earlier - which is partially resolved by 
evolving h in a particular way when using MLS approximation - we period- 
ically update the particle volume predicted by (12) with J <j>i dx if there is 
significant difference between these two quantities. To be more specific, the 
particle volume Vi is updated if \Vi — J <pidx\/Vi > l.Ox 10 ~ 3 . 

We first ran a simulation with linearly complete MLS shape functions. The 
underlying univariate function, w, was selected to be a Wendland function 
with C 4 -smoothness. The smoothing-length was evolved by taking a time 
derivative of the relationship hi oc Vi and integrating it alongside the other 
equations, the constant of proportionality was chosen to be 2.0. The result 
is shown in Fig. 2. The agreement with the analytical solution (solid line) 
is excellent, especially around the contact discontinuity and the head of the 
rarefaction wave. Next, we constructed RBF shape functions. As we mentioned 
in Sect. 3, for this one-dimensional problem we employ cubic B-spline because 
they constitute discrete bi-Laplacians of the shifts of the globally supported 
basis function, ip = | • | 3 . The result of this simulation is shown in Fig. 3. 
Again, the agreement with the analytical solution is excellent. 

In the introduction an SPH simulation of the shock tube was displayed 
(Fig. 1). There, we integrated (4)-(6) and h was updated by taking a time 
derivative of the relationship hi = 2.Qrm/pi. To keep the comparison fair, 
the same initial condition, particle setup and dissipative terms were used. As 
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Fig. 2. Shock tube simulation (t=0.2) using linearly complete MLS shape functions. 
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Fig. 3. Shock tube simulation (t=0.2) using cubic B-spline shape functions. 



previously noted, this formulation of SPH performs poorly on this problem, 
especially around the contact discontinuity. Furthermore, we find that this 
formulation of SPH does not converge in the Loo-norm for this problem. At 
a fixed time (t = 0.2), plotting number of particles, N, versus Loo-error in 
pressure, in the region of the computational domain where the solution is 
smooth reveals an approximation order of around 2/3, attributed to the low 
regularity of the analytical solution, for the MLS and RBF methods, whereas 
our SPH simulation shows no convergence. This is not to say that SPH can not 
perform well on this problem. Indeed, Price [11] shows that, for a formulation 
of SPH where density is calculated via summation and variable smoothing- 
length terms correctly incorporated, the simulation does exhibit convergence 
in pressure. The SPH formulation we have used is fair for comparison with 
the MLS and RBF methods since they all share a common derivation. In 
particular, we are integrating the continuity equation in each case. 
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To conclude, we have proposed a new set of discrete conservative variation- 
ally-consistent hydrodynamic equations based on a partition of unity. These 
equations, when actualised with MLS and RBF shape functions, outperform 
the SPH method on the shock tube problem. Further experimentation and 
numerical analysis of the new methods is a goal for future work. 
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